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Stochastic models of reaction-diffusion systems are important for the study of bio- 
chemical reaction networks where species are present in low copy numbers or if reac- 
tions are highly diffusion limited. In living cells many such systems include reactions 
and transport on one- dimensional structures, such as DNA and microtubules. The 
cytoskeleton is a dynamic structure where individual fibers move, grow and shrink. 
In this paper we present a simulation algorithm that combines single molecule simula- 
tions in three-dimensional space with single molecule simulations on one-dimensional 
structures of arbitrary shape. Molecules diffuse and react with each other in space, 
they associate to and dissociate from one-dimensional structures as well as diffuse 
and react with each other on the one-dimensional structure. A general curve embed- 
ded in space can be approximated by a piecewise linear curve to arbitrary accuracy. 
The resulting algorithm is hence very flexible. Molecules bound to a curve can move 
by pure diffusion or via active transport, and the curve can move in space as well as 
grow and shrink. The flexibility and accuracy of the algorithm is demonstrated in 
four numerical examples. 
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I. INTRODUCTION 



It has been established that modeling biochemical reaction networks with a macroscopic 
deterministic model, for instance with ordinary differential equations (ODEs), can be unsuit- 
able when some species are present in low copy numbers or if reactions are highly diffusion 
limited^. Such systems may only be accurately modeled by a discrete stochastic model. 
Usually two levels of stochastic modeling are considered: the mesocopic level and the mi- 
croscopic level. 

The mesoscopic scale is described by the chemical master equation (CME) which is a 
discrete, stochastic model. Molecules are treated as discrete entities but assumed to be 
well stirred with a uniform distribution within the simulation volume. The copy number 
of each species is computed but the position of individual molecules is not tracked. Exact 
trajectories of such a system can be generated with the stochastic simulation algorithm 
(SSA) by Gillespie^. 

Systems for which the condition of spatial homogeneity does not hold can be modeled at 
the mesoscopic level by subdividing the domain into smaller compartments, where in each 
compartment the well stirred condition is assumed to hold. Diffusion can now be modeled 
as discrete jumps between the compartments. The governing equation of the system is the 
reaction-diffusion master equation (RDME) and exact trajectories can be computed with 
the next subvolume method (NSM) 2 , which is an efficient implementation of the SSA in a 
spatial setting. Several software packages have been developed for simulating systems at 
the mesoscopic level, including MesoRD 31 , URDME 3 and STEPS 21 . The spatial accuracy 
of the RDME depends on the size of the voxels. If the voxels are too large we might not 
resolve the spatial dynamics of the system accurately, but on the other hand it has also been 
shown that when the size of the voxels becomes too small the accuracy deteriorates^. By 
allowing the reaction rates in the system to depend on the size of the voxels in a general 
way, accuracy can be retained, but at most up to a critical lower bound of the voxel size, 
at which point the model breaks down regardless of how the reaction rates are choserP^. 
The mesoscopic model may therefore be unsuitable for systems where a very high spatial 
accuracy is required to capture important dynamics. 

The microscopic level, or microscale, usually refers to models where the positions of single 
molecules are tracked. Either we compute the continous position in space or let the individual 
molecules move on a lattice, as in Ref. [TTJ Reactions occur with some probability per unit 
time when molecules collide or when they are within some prescribed reaction radius. There 
are two models that have been used extensively to model systems at the microscale: the 
Smoluchowski model 1? and the Doi modeP^. In this paper we will only be concerned with 
the Smoluchowski model. 

In Ref. US] they show that fine-grained spatial correlations can affect the macroscopic 
behaviour of a system. Furthermore, the dynamics of the system cannot be captured on the 
macroscopic or mesoscopic level but rather require a microscopic level of simulation. This 
is due to highly diffusion-limited reactions and fast rebindings that the mesoscopic model 
cannot resolve without hitting the critical lower bound of the voxel size. This motivates 
development of efficient methods for simulating systems at the microscopic level. 

The Green's function reaction dynamics (GFRD) algorithm^ is an efficient method to 
generate accurate trajectories of the Smoluchowski model. A full iV-body problem is divided 
into many 1- and 2-body problems that can be simulated independently by choosing the time 
step in the algorithm sufficiently small. For these subproblems the Smoluchowski equation 
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can be solved analytically or numerically. Improvements of GFRD have been developed 
in e.g. Refs. [15] and [T71 Sampling random numbers from the analytical solution of the 
Smoluchowski equation is complicated and expensive, but in Ref. [T7]this step is simplified 
by solving the Smoluchowski equation using first order or second order operator splitting 
schemes. In Ref. [T5| the algorithm is made exact in both space and time by combining 
GFRD with the first-passage kinetic Monte Carlo (FPKMC) algorithm developed in Ref. 
[T8l and this improved algorithm is called eGFRD. The main idea in eGFRD is to introduce 
protective spheres around single molecules and pairs of molecules and then sample exact 
times for the molecules to exit the spheres. 

Software that implements eGFRD is eCelP^. Other software packages that simulate 
systems with the Smoluchowski model as the target model are SmoldyrP^, MCelP^ and 
Spatiocyte^. Both Smoldyn and MCell discretize time while GFRD is exact in time, thus 
making them less accurate. On the other hand, Smoldyn and MCell can be significantly 
more efficient for larger systems. In Spatiocyte molecules move on a lattice of densely packed 
spheres. 

Instead of simulating the system with a purely microscopic method one can combine 
the microscopic level and the mesoscopic level in a hybrid methocP^HSl gy restricting the 
microscale simulations to the species or the parts of space that are necessary for capturing 
the important dynamics of the system, computational time can be saved without losing a 
significant amount of accuracy. 

In addition to the observed spatial heterogeneity in systems where molecules move only 
by pure diffusion in free space, there can be spatial effects in systems due to interactions with 
internal lower dimensional structures in the cells. This could for instance be binding to and 
reactions on curved surfaces or one-dimensional structures such as DNA or microtubules. In 
Ref. 122] a method for simulating interactions with curved surfaces at the microscopic level is 
developed, and in Ref. [25] they develop an FPKMC method for simulating active transport 
due to chemical potentials on embedded lines of fixed length and with fixed configuration. A 
mesoscopic method that simulates interactions between molecules in free space and lines is 
considered in Ref. [26] and a mesoscopic method incorporating active transport is suggested 
in Ref. 123 

In this paper we have developed an algorithm that combines single molecule simulations 
in three-dimensional space with single molecule simulations on embedded curves. Molecules 
in some general three-dimensional volume, for instance a sphere or a cylinder, are simulated 
with the GFRD algorithm. They can react with the embedded curves according to the 
Smoluchowski model. Molecules attached to a curve can also diffuse and react with each 
other according to the one-dimensional Smoluchowski equation. In addition to pure diffusion, 
we simulate active transport along the curves by letting the molecules bound to the curves 
move deterministically in one direction. The curves can have an arbitrary shape. This 
flexibility is obtained by approximating a continuous curve with non-constant curvature by 
straight linear segments, and this approximation can be made arbitrarily accurate by making 
the discretization of the curve finer and finer. It is also known that many one-dimensional 
structures in cells are highly dynamic^. A microtubule may for instance both move in space 
as well as grow and shrink. We take those dynamics into account by letting the embedded 
curves move according to some general random process as well as letting them grow and 
shrink over time. 

The following sections are organized as follows. In Sec. [IT] we review the GFRD method 
in one dimension as well as three dimensions. In Sec. |III| we describe the algorithm for 
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coupling the one- dimensional simulations with the three-dimensional simulations, as well as 
simulation of moving lines embedded in space. In Sec. IV we demonstrate the flexibility 
of the algorithm with four examples. In the first example we verify the accuracy of the 
method by comparing with mesoscopic simulations and theoretical predictions. In the second 
example we consider two spirals embedded in a sphere. In the third example we consider 
several straight lines embedded in a spherical domain. Molecules can bind to the lines and 
are then moved by active transport along the lines. The lines move in space according to 
a random process, and we show how the active transport affects the distribution of the 
molecules in space. In the final example we consider lines that grow and shrink. We show 
how the number of molecules attached to a line will vary with the total length of the lines. 

All parameters are given in Si-units. 



II. SIMULATIONS IN THREE AND ONE DIMENSIONS 

We consider a system of N molecules, where the molecules can undergo bimolecular 
or unimolecular reactions. The evolution of the system is modeled by the Smoluchowski 
equation with a mixed boundary condition at the reaction radius of the molecules. In 
principle one could write down the Smoluchowski equation for the full system, but due to 
the high dimension it would not be tractable to solve. In GFRD^the main idea is to divide 
the full A-body problem into smaller subproblems consisting of the simulation of single 
molecules and pairs of molecules. The subdivision is done in two steps. First we find the 
nearest neighbor of every molecule. Molecules that have each other as nearest neighbors 
form pairs and the remaining molecules will be updated as single molecules. Secondly we 
choose a time step At such that no molecule in a pair is likely to react with any other 
molecule except its nearest neighbor and molecules that are not in a pair should be unlikely 
to react with any other molecule. The full A-body problem can now be approximated by 
smaller, independent, subproblems during the time step At. 



A. Updating positions of single and pairs of molecules 

The positions of single molecules are updated by sampling from a normal distribution 
with standard deviation \/2DAt in each direction, where D is the diffusion constant of 
the molecule. A pair of molecules, where the molecules have the positions xi„ and X2 n at 
time t n , are updated according to the PDF p(xi, x 2 , t|xi„, x 2 „, t n ) obtained by solving the 
Smoluchowski equation for a pair of molecules, given by: 

Pt = D 1 A Xl p + D 2 A X9 p, (1) 

where Di and Xj are the diffusion coefficient and position of particle i. By making a change 
of variables as suggested in Ref. [16] this equation can be split into two equations: one for 
the relative position of the molecules and one for a weighted mean of the positions 

Y = v^X! + v/£)i/£) 2 x 2 , y = x 2 - Xl . (2) 
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Now p can be written as p(xi, x 2 , t|x in , x 2n , t„) = pv(Y, t|Y„, t„)p y (y, t\y n , t n ), and the 
equations become 

d t p Y = DAyPy, d t p y = DA y p y , 

where D = D1+D2 is the sum of the diffusion coefficients of the molecules. The motion in the 
Y direction is in free space and will be described by a normal distribution, but the equation 
for the y coordinate will have an additional boundary condition at the reaction radius of 
the molecules to model associations between molecules. By considering the equation in a 
spherical coordinate system, r = (r,9,(f>), the boundary condition becomes: 



47RX 2 L> 



dp r 

Or 



(3) 



where a is the reaction radius of the molecules, r = ||r|| and k r is the intrinsic association 
rate. The initial condition is given by 

p r (r,t\r n ,t n ) =S(t- r n ) 

and p T (r,t\r n ,t n ) vanishes as r — > 00. 

This equation can be solved analyticall y ^ 16 * 29 !, but due to the complexity of the analytical 
solution we use the method described in Ref. [T71 There the full equation is solved in two 
steps using a first order splitting scheme (or if higher accuracy is needed with a second order 
splitting scheme). This simplifies the sampling of new positions for the molecules, at the 
price of an additional error. This error is shown to be small in Ref. [TTJ The first step of 
the operator splitting scheme is to solve for the radial part of the equation (in a spherical 
coordinate system rotated such that r n = (r n , 0, 0)), given by: 



9 tP = D f— + (4) 
r \ d 2 r r dr J 

with the initial condition p r (r, t n \r n , t n ) = 8(r — r n ) and with the boundary condition 
given by Eq. (J3]). The analytical solution is derived in Ref. [3Ql The second step is to solve 
for the angular part: 



D ( 1 d ( . dp e \ 1 d 2 p e \ 

with initial condition given byp e (#,t„|#„,0 n ,t n ) = 5(9)/r 2 sm(9). The solution to Eq. Q 
will be independent of 0, and the distribution of <fi will therefore be uniform on the interval 
[0, 27r). It is possible to solve Eq. ^ analytically and the solution is given in the form 
of an infinite series, where the terms are products of Legendre functions and exponentials. 
Efficient strategies to compute the analytical solutions of Eqs. Q and ^ are discussed in 
Ref. [3TJ An efficient method to generate random numbers from the PDF pe is suggested in 
Ref. |32J 



5 



Note that the probability for two molecules to survive until time t n + At is given by 

" tM dp r 



S(tn + At\r n , t n ) = l- A7ia 2 D j 



dr 



dr. 



From S we can sample the exact time when the pair of molecules react (cf. Ref. [TB]) . Thus, 
the first step is to check whether the two molecules react before t n + At. If they do, we 
execute the reaction at the sampled time, and if they do not react we sample new positions 
for each molecule by 

1. Sample r n+ i from the PDF p r . 

2. Sample # n+ i from the PDF p e with r = r n+1 . 

3. Sample n+ i from a uniform distribution on the interval [0, 2ir). 

4. Transform r n+ i back to y n +i- Sample Y n +i from a normal distribution and then solve 
Eq. §2§ to obtain x 1(n+ i) and x 2(n+ i). 

Also molecules in one dimension can react. For pairs of molecules we solve the one- 
dimensional Smoluchowski equation 

d tPs = (6) 
with the boundary condition at the reaction radius 

D-^-\ s=tJ = k r p s (s = a,t\s n ,t n ). (7) 

to obtain the PDF p(s, t\s n , t n ) for the distance s at time t between the two molecules, given 
that the distance was s n at time t n . The analytical solution to Eq. ^ with boundary 
condition given by Eq. ^ is derived in Ref. [2H1 

Both molecules in one dimension and molecules in three dimensions can dissociate. Con- 
sider a reaction A — > B + C. The time until an A-molecule dissociates is sampled from 
an exponential distribution with mean kd- After a molecule A has dissociated the B- and 
C-molecules are placed at a distance equal to the reaction radius. 



B. Molecules close to a curved boundary 

Most biologically realistic geometries are fairly complex. Cells can be close to spherical 
or have a shape reminding of a spheroid. Eukaryotes typically have complex internal two- 
dimensional structures. One flexible and efficient way to deal with such non-regular surfaces 
is to discretize the surface with an unstructured mesh. 

In this paper we do not consider volumes with internal two-dimensional structures, and 
are consequently only concerned with the outer boundaries of the volume. The outer bound- 
ary is approximated by a triangular surface mesh. The vertex of each triangle is also a center 
vertex in the dual mesh generated from the primal mesh. Consider a center vertex ipi. This 
center vertex is the vertex of several triangles T\, . . . , Tm in the mesh, and to each dual ele- 
ment we associate one linear plane defined by the center vertex ipi and with normal given by 
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the average of all the normals nj, . . . , njj to the triangles 7\, . . . , Tm- We have now assigned 
one linear plane to each dual element in the mesh, and these planes are an approximation 
of the surface. 

Now, assume that a molecule is close to one of the planes approximating the surface. 
Choose a time step At small enough such that the interaction between the molecule and 
the plane can be approximated by the interaction between an infinite flat plane and a 
molecule during At. Then consider the change of coordinate system where the molecule 
diffuses parallel to the plane in two directions and along the normal of the plane in the third 
direction. The molecule can then be updated according to normal diffusion in the first two 
directions and in the third direction by sampling from a normal distribution and then, in 
case the molecule ends up on the wrong side of the plane, simply reflecting that coordinate 
in the plane. 

This approach gives great flexibility to the algorithm and enables us to simulate molecules 
in different complicated geometries. The details of this approach is described in Ref. [221 
Another method to simulate molecules close to an absorbing or reflecting boundary has been 
suggested in Ref. |33l 

Note that by making a locally planar approximation we can also include reactions between 
molecules and surfaces by solving the one-dimensional Smoluchowski equation for molecules 
that are close to the surface, as done in Ref. |22| In this paper, however, we only consider 
reflecting outer boundaries, but the methodology for simulating interactions with reactive 
surfaces (outer boundaries as well as internal surfaces) developed in Ref. [22] can be combined 
with the method described in this paper. Furthermore, we could also combine this method 
with the hybrid method coupling microscopic and mesoscopic simulation algorithms, also 
developed in Ref. 

III. EMBEDDED ONE-DIMENSIONAL MANIFOLDS 

We consider a continuous and differentiable curve F embedded in space, parametrized by 

r(s) = (a(s), b(s), c(s)) , < s < 1. 

T can be approximated by a piecewise linear curve 7. Denote the linear segments of 7 by 
7i,...,7tv- The segment 7, is defined by the points pj and pf. The points are chosen such 
that p\ = pl_i and thus all segments are connected. The curve V can be approximated to 
arbitrary accuracy by letting N — > 00 in such a way that \\pj — p?|| — > 0. 

We want to simulate reactions between molecules in space and the curve as well as 
diffusion and reactions between molecules on the curve. Furthermore, we want to simulate 
moving curves that grow and shrink. 

A. Interactions with the curve 

If a molecule A is close to 7, we will choose a time step At such that the molecule is 
unlikely to interact with more than one of the linear segments during that time step, and 
such that the linear segment can be approximated, to high accuracy, by an infinitely long 
straight line during At. The problem of simulating a molecule close to an infinitely long 
straight line can be reduced to a two-dimensional problem, which can be seen by considering 
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the position of the molecule at time t n , x n , in cylindrical coordinates r = (r, 9, z), with the 
^-coordinate in the direction of the line. The diffusion of the molecule in the ^-coordinate is 
now independent of the diffusion in the r and 9 directions and the interaction with the line. 

Assume that the closest linear segment to the molecule is 7$. Now consider the sphere S 
with center at x n and radius R = min{||x n — p}\\, ||x n — p 2 ||}. By choosing the time step At 
such that the molecule A is unlikely to diffuse outside of the sphere S during the time step, 
the requirement that the molecule is unlikely to react with more than one linear segment 
is satisfied. If the diffusion coefficient of A is Da, the molecule will, on average, diffuse the 
distance \/6D^At during a time step At. Thus, we choose the time step to be 

R 2 

At 



K ■ 6D A 



where K is a constant chosen sufficiently large in order to ensure that the probability for the 
molecule to diffuse outside of the sphere is small. Note that if the molecule is really close to 
pj or p 2 the time step could become very small, making the simulation too slow. To avoid 
this problem we have to choose a smallest value that At can take, thus introducing an error. 
However, with a carefully chosen smallest value of At that error will be small. 

The interaction between the molecule and the linear segment can now be approximated by 
the two-dimensional Smoluchowski equation during At. The two-dimensional Smoluchowski 
equation for the relative position (in cylindrical coordinates) is given by: 

Pt = DAp, (8) 

with the boundary condition 

dp r 



2ttctD^- 
or 



k r p r (r = a,t\r n ,t n ). (9) 



Again we solve the equation in two steps with a first order operator splitting scheme. 
The radial part is given by: 

»*r=(^ + ^) (10) 



dr 2 r dr t 

with the boundary condition given by Eq. pi), the initial condition given by p r (r, t n \r n , t 



5(r — r n ) and with p r vanishing at infinity. The solution of Eq. (10) with the boundary 
condition given by Eq. ^ is solved analytically in Ref. |29l The analytical solution is 



1 



00 



p r (r,t\r ,t n ) — — / exp(—Du(t — t n ))C(u,r,k,k r )C(u,ro,k,k r )vdu, 

Z7T 







where k = 2iraD and the function C is defined by 

Jo(ur)(kuYi(au) + hY Q (au)) — Yo(ur)(kuJi(au) + hJo(au)) 



C(u, r, k, h) 



{{kuY 1 (au) + hY (au)) 2 + (kuJ^au) + hJ {au)) 2 )z 



Here Jo and J\ are Bessel functions of the first kind and Yq and Y\ are Bessel functions 
of the second kind. Despite having the analytical solution available we choose to compute 
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the solution using a finite difference method, as the analytical solution is fairly complex 
and expensive to evaluate. The solution could also be tabulated as described in Ref. I3T1 
Tabulating the solution is more efficient but has the drawback of introducing an additional 
error as well as consuming more memory. 
The second step is to solve for 9 and z 

a n ( 1 d 2 p 6z d 2 p 0z \ 

dpe z = D + , (11) 



with initial condition 



r 2 d9 2 dz 2 



(a + \ 6 ( 9 ) 6 ( z ~ z r> 
p dz (6,z,t n ) = 



Again, p vanishes at infinity. The solution of Eq. (11) is a Gaussian in both the ^-direction 
and the ^-direction. 

Molecules can also dissociate from the lines. The time until a dissociation is assumed 
to be exponentially distributed with mean equal to the intrinsic dissociation rate, and after 
a dissociation from a line the molecule is placed at a distance equal to the reaction radius 
between the molecule and the line. As we are using a finite difference scheme to compute 



the solution of Eq. ( 10 ) we are in fact forced to place the dissociating molecule at a slightly 
larger distance from the line than the reaction radius, in order to have discretization points 
between the molecule and the line. This introduces a small error, but it can be made 
arbitrarily small by placing the molecule closer and closer to the line while increasing the 
number of discretization points used in the finite difference scheme. 



B. Diffusion and reactions on moving curves 

A new position at time t + At, x n+1 , of a single molecule moving by pure diffusion, bound 
to T at position x n at time t, where x n e 7^ is computed by sampling a distance d from a 
one-dimensional normal distribution with standard deviation \/2DAt and then updating 

x n+l = x n + ^TT~! ^~TTr (12) 

IIP? - Pi II 

It is possible that x n+ i is no longer on the segment If this is the case we project x n+ i to 
the closest point on 7. The time step At should thus be chosen such that \j2DAt is small 
in comparison to ||p 2 — pf||, so that the molecules do not diffuse a long distance compared 
to the distance of the segments. Since we are projecting the molecule onto the line we are 
introducing a small error and the effective diffusion will be slightly slower than the exact 
diffusion process. However, if the angles between the linear segments are small, this error 
will also be small, and it is therefore important to ensure that the discretization of V is fine 



enough to make the angles small enough. Note that the process defined by Eq. (12) can 
be modified to account for other types of motion on curves, such as active transport. For 
instance, d could be sampled from another distribution than the normal distribution, or be 
completely deterministic. Molecules on the same line segment can react according to the 
PDF obtained by solving Eq. ^ with boundary condition given by Eq. ([7]). 

Many one-dimensional structures in cells are not stationary, but rather move around in 
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space. They can also grow and shrink. We simulate this process with an operator splitting 
scheme. Choose a time step At sp u t . First we keep the lines constant in space and propagate 
the molecules according to the algorithm described above for the full time step A£ S piit- The 
second step is to move the lines according to some transformation T(p,t) : M. 3 x M + — > M. 3 
in space, for the full time step At split . T can both move the line in space but also change 
the length of the line. We apply T to all the points, such that the linear segments at time 
t + Aigpiit are defined by the pair of points 

(T(p{, Atsput), T(pl, At split )), . . . , (T(pjv, At split ), T(pjf, At split )). 

The position of the molecules are now no longer guaranteed to be on the lines. Each 
molecule is therefore projected down to the closest point on its respective line. For high 
accuracy we should therefore choose the time step At sp i it such that each molecule is projected 
only a short distance. After each step it is possible that some molecules on the line that 
were previously separated by some small distance are now overlapping. If this is the case 
they are moved apart and placed at a distance equal to their reaction radius. Molecules in 
space that overlap with the lines are also moved such that they are at a distance from the 
line equal to their reaction radius. We then repeat this process until the final time tf. If 
higher accuracy is needed it would be possible to consider higher order splitting schemes, 
such as Strang splitting^. 

The algorithm outlined above is summarized in Algorithm 1. Note that this is only a 
summarized version of the GFRD algorithm. There are some special cases that require 
extra care, such as pairs of molecules that are also close to a boundary. For a more detailed 
description of how to handle such cases the reader is referred to Refs. ITST - TITI 

IV. NUMERICAL RESULTS 

We show the practical use of the algorithm by studying both examples with straight 
lines through a three-dimensional domain as well as lines with non-constant curvature. We 
demonstrate how to use the algorithm to simulate active transport by letting molecules 
attached to a line move with a deterministic speed in a predetermined direction. Finally, 
we let the lines move in space as well as grow and shrink. 

A. Polymer with road blocks embedded in a cylinder 

We consider a cylinder with radius R = 10~ 6 , height H = 2 ■ 10~ 6 and aligned along the 
x-axis. Inside the cylinder we consider a line of length H, for instance modeling a polymer, 
going through the center of the cylinder. The polymer has a reaction radius 77 = 10~ 9 . 
Molecules of species A can associate to the polymer to form molecules A cy \. Molecules A cy \ 
can dissociate from the polymer to become species A. Thus, we consider the reaction 

A — - A cy \. 

kd 

where A is a molecule diffusing in three-dimensional space and A cy \ is a molecule bound to 
the polymer. The diffusion constant of A is Da = 1CT 12 and the diffusion constant of A CJ \ 
is Da = 10~ 14 . The intrinsic reaction rates are chosen to be k r = 10 -11 and = 50. The 
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ALGORITHM 1: GFRD with moving one dimensional structures 



1: Initialize the system at the time t = and choose a final time tf. Choose a time step 

A4piit for the operator splitting. 
2: Divide the set of all molecules into subsets of one or two molecules such that pairs are 

formed by molecules that are each others nearest neighbor. Set ti oc = 0. 
3: while ti oc < AtspHt do 

4: Choose a time step At < At sp ii t — t\ oc such that molecules are unlikely to interact with 
more than one other molecule during that time step. At should also be chosen such that 
no molecule is likely to interact with a two-dimensional boundary or a one-dimensional 
curve and a neighbor during the time step. To fulfill this requirement certain pairs and 
single molecules may have to be updated in several smaller steps during At. 
5: Pairs of molecules are updated according to the PDF obtained by solving the three- 
dimensional Smoluchowski equation ([I]), as described in Sec. II A 
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6: Molecules close to a line are updated according to the PDF obtained by solving the 
two-dimensional Smoluchowski equation ( |To] ), as described in Sec. Ill A 



7: Molecules on a line are updated according to the PDF obtained by solving the one- 



dimensional Smoluchowski equation (|6j), as described in Sec. II A 

Molecules that dissociate from the cylinder at time r < At are placed at the reaction 
radius of the molecule and the cylinder. If a molecule in free space dissociates we place 
the products at a distance equal to the sum of the reaction radii of the two molecules. 

Update t\ oc = t\ oc + At. 
end while 

Update the linear segments according to the transformation T. 
Update the time t = t + At sp ii t . 
Repeat (2)-(12) until the final time tf. 



problem setup is shown in Fig. [TJ 

Due to the cylindrical symmetry of the problem, the reactions with the polymer can be 
reduced to a two-dimensional problem. The average time Tbmd until a molecule A reacts with 
the polymer can thus be computed using a formula for association times on disks, derived 
in Ref. HOI We can therefore compare the microscopic simulation with a corresponding 
mesoscopic simulation of the reaction 

j, meso 

^4 \ ^-cyl 

Ltncso 

K d 

where fc™ eso = 1/Tbind and, due to the detailed balance condition, fc™ eso = kdk™ eso / (AdiskK) . 
Here Adisk = nR 2 is the area of the face of the cylinder. In Fig. [2] we show that for 
this simple problem, the microscopic simulations agree well with spatially homogeneous 
mesoscopic simulations, performed with the SSA. 

Next we consider a more involved example. The geometrical setup is the same as above, 
but now we let the molecules of species A cy \ react with an operator site at the center of the 
line, thus considering the line a coarse model of DNA. The operator site is denoted by B 
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FIG. 1: The problem setup. Molecules of species A (blue dots) diffuse in the cylindrical volume. 
They react with the line and transforms into the one-dimensional species A cy \ (red dots). The 
line (green) extends throughout the domain through the center of the cylinder. The outer 
boundary is approximated by piecewise linear planes and is reflective. 
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FIG. 2: The average number of A cy \ at the microscopic level and the mesoscopic level. The 
mesoscopic simulation agree well with the microscopic simulation. The average is taken over 20 
trajectories and we start the simulation with 250 A-molecules uniformly distributed in the 

domain. 



and we consider the reactions 



A A cy \ 

kd 

A cy i + B — > B a 
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FIG. 3: The average time until an ^4-molecule binds to the operator site decreases with the 
distance between the road blocks. As the distance becomes greater than the average sliding 
distance of an A cy \, the decrease levels off. The parameters in the simulations are k r = 10~ n , 
kd = 50 and D = 10~ 12 for all species. The simulation is started with 250 A-molecules in free 

space. 



We assume that the operator site is absorbing and of radius 10~ 9 . On both sides, at the 
distance Z r b from the operator site, we place a road block. Road blocks are modeled by 
stationary, non-reactive (but reflective) molecules of radius 10~ 9 . Thus, molecules binding 
to the DNA at a distance larger than from the operator site, will not be able to find it since 
molecules cannot diffuse through the road blocks. In Ref. [26] a similar problem is simulated 
with a mesoscopic method. The microscopic parameters cannot easily be transformed to 
corresponding mesoscopic parameters, but the qualitative behaviour can be compared. 

We compute the average time, T act i V e, until an A cy i-molecule binds to the operator site, 
and this time will depend on the length Z r b between the road blocks. In Fig. [3] we show 
that this time decreases with increasing length between the road blocks, as was shown in 
Refs. [26] and [35j The decrease in the binding time levels off as / r b becomes greater than the 
average distance that a molecule slides on the DNA before dissociating. 



B. Spirals embedded in a sphere 

One-dimensional structures in living cells are often not well approximated by straight 
lines due to their curvature, and they may also have a complicated spatial distribution 
relative to each other. It is therefore important to be able to simulate molecules reacting 
with and on lines of non-constant curvature. To demonstrate the flexibility of the method 
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FIG. 4: We have simulated association to and dissociation from two spirals (yellow and red lines) 
embedded in a sphere. Molecules in space diffuse and react with the spirals, and molecules on the 
spirals diffuse and react with each other. Note that the molecules that are free in space are not 

plotted in the figure above. 

we therefore consider the following system of reactions 

A ^ A cyh B B cy \ 

^■Cyl + B C y\ V— 1 C C yl, 

k d 

where A cy \, B cy \ and C cy \ are molecules bound to one of the lines and A and B are molecules 
diffusing in three-dimensional space that can react with the lines. We will consider two 
spirals I\ and T2 along the y-axis, parameterized by 

T 1 {a{s),b{s),c{s)) = (-3 • 1(T 7 ,0,0) + (2.5 • 10~ 8 s, r c cos(2vrs), r c sin(2vrs)) 
r a (o(s), b(s), c(s)) = (-3 • 1(T 7 , 0, 0) + (2.5 • 10' 8 s, -r c cos(2vrs), -r c sin(27rs)), 

with < s < 3 and r c = 3 • 10~ 7 . Each spiral is divided into 30 linear segments. The setup 
is shown in Fig. |4j 

We start with 150 molecules of species A and 150 molecules of species B, with initial 
positions sampled from a uniform distribution. The reactions rates are chosen to be k\ = 

= 10~ n , k\ = k\ = 50, k^ = 10~ 6 and k\ = 20. The diffusion rate for the molecules in 
space is 10 -12 and the diffusion rate for the molecules on the spirals is 10 -14 . We simulate 
the system for 2 seconds and the evolution of the system is shown in Fig. [5| 

C. Active transport on moving microtubules 

Molecules in living cells do not only move by normal diffusion but can also be actively 
transported. For instance, intracellular transport on the cytoskeleton, a structure in the cell 
composed of fibers such as actin and microtubules, is an important mechanism for directed 
transport of cargo in the celP^. Cargo attach to the fibers and are then transported by 
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FIG. 5: The time evolution of the total number of molecules of the different species in the system. 



molecular motors: myosins along actin and dynein and kinesins along microtubules. 

We will consider a much simplified model of active transport from the cytosol towards the 
center of a cell. Molecules bind to microtubules, modeled by straight lines passing through 
the origin. We assume that the transport is of constant speed in the direction towards the 
origin and we do not explicitly model the molecular motors. Other more complex and more 
biologically accurate models of active transport could be incorporated into the method by 



modifying Eq. (12) in a suitable way. The model used in this example is not meant to be 



a very accurate model of such processes, but is rather meant to serve as an example of the 
flexibility of the algorithm. 

We initialize the system by sampling 25 points, pi, . . . , P25, uniformly on a sphere of 
radius R = 1CT 6 . 25 straight line segments through the origin are then defined by the 
pairs (pj, — Pi), i = 1, ... ,25. We start with 100 A-molecules in space, with initial positions 
sampled from a uniform distribution. The A-molecules can then associate and dissociate 

to the microtubules through the reaction A A cy \. A- molecules in space diffuse with the 

diffusion constant D = 10~ 12 , but molecules on the microtubules move deterministically 
towards the center of the sphere with constant speed \/2 ■ 10 _14 At/40. The microtubules 
are rotated A9 and A<p radians around the x- and y-axis per time step At, where 



A9 = 27rAtXi 
A0 = 2vrAtX 2 . 



(13) 



Xi and X 2 are random numbers sampled from a normal distribution with standard deviation 
1. At S piit is chosen to be 0.01. Thus, the microtubule U defined by the points (pi, — p«) at 
time t, will be updated at time t+ At by rotating pi and — pi around the x- and y-axis AO and 
A(f> degrees, with AO and A<p given by Eq. (13). The trajectory of one of the microtubules 
is plotted in Fig. [6] and two snapshots of a simulation, one at t\ = 0.05s and one at the final 
time tf = 1.0 s are shown in Fig. [7j 

In Fig. [8] the distribution of molecules in the radial direction is plotted. The initial 
distribution is uniform in space. Due to the transport towards the center, this distribution 
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FIG. 6: The trajectory of one of the microtubules. The position of the x-, y- and z-coordinates of 
one of the endpoints of the microtubule is plotted. The microtubule is rotated according to Eq. 



(13) 





(a)t = 0.05s 



(b)i = 1.0s 



FIG. 7: Lines are yellow and the blue dots are molecules bound to the lines. Molecules in free 
space are not plotted. Molecules on the lines move by active transport towards the center of the 



sphere and the lines move in space according to Eq. (13) 



is shifted and more molecules are close to the center of the sphere at the final time tf = 1.0s. 
The reaction rates in the simulation are k r = 10 -11 and kd = 10 3 . 



D. Growing and shrinking microtubules 



Microtubules are highly dynamic and can, in addition to moving around in space, both 
grow and shrink. Therefore, as a final example, we will consider fibers embedded in a sphere 
and whose lengths can both increase and decrease. 

The sphere has a radius R = 10~ 6 , with two straight fibers passing through the origin. 



The fibers are initially of length I = 2R. The lines rotate in the same way as in Sec. IV C 
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FIG. 8: The radial distribution of the positions of the molecules. The molecules have a uniform 

initial distribution in the sphere (dashed line). The distribution at the final time tt = 1.0s is 
shifted towards the center of the sphere. The distribution at the final time, plotted in the figure 

above, is based on the mean of three trajectories. 



but will also grow and shrink. After each time step, the length of the line is updated, and 

molecules are projected down to the closest point on the line. Note that the closest point 

could in fact be the endpoint of the line if the molecule was close to the end of the line 

before the length was updated. One possibility is to let such molecules be transformed to 

molecules in free space, but in this example they will stay bound to the line. Now, assume 

that the length of the line at time t n is given by l(t n ). Then the length at t n + At is given 
by ^ • ■ 

l(t n + At) = l(t n ) + X, 

where X is a random number sampled from a normal distribution with mean and standard 
deviation y/2DiAt, with the restrictions that l(t n + At) < 2R and l(t n + At) > 1CT 8 . Thus, 
if the line is denned by the points (p n , — p n ) at time t n , it will be defined by the points 
l(t n + At)//(t„) • (p n , -p n ) at time t n+1 . 

Intuitively one would expect that when the length of the fibers decrease, the number of 
bound molecules will also decrease, under the condition that the fibers grow and shrink at 
a time scale that is longer than or similar to the time scale for the system to reach steady 
state. Thus, if we choose parameters such that this condition is fulfilled, we would expect 
to see a correlation between the number of molecules of species A cy \ and the total length of 
the two lines. 

We start the simulation with 600 molecules of species A uniformly distributed in space. 
The diffusion constants are Da = 10~ 12 and D^ cyl = 10~ 14 . In Fig. [9] we plot the total 
length of the two lines and the total number of molecules bound to a line. In Fig. [9] a 
trajectory of the system, with Di = 10~ 12 , is plotted together with the sum of the length of 
the lines. As expected the two trajectories are clearly correlated. 
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FIG. 9: The number of ^4 cy i-molecules is plotted together with the sum of the length of the two 
lines. The two time series are correlated with correlation coefficient 0.713 and a 95% confidence 
interval of the correlation coefficient is given by [0.669,0.752]. 



V. CONCLUSIONS 

We have described an algorithm to simulate biochemical systems in complex geometries 
with embedded and dynamic one-dimensional manifolds. Molecules in space can associate to 
and dissociate from the embedded curves. Molecules in space and on the lines are simulated 
using the GFRD algorithm. 

The accuracy of the algorithm depends on the discretization of the lines and the time step 
chosen in the algorithm. The simulation of the molecules and the spatial transformation of 
the lines are split in time according to a first order operator splitting scheme and the error 
will therefore depend on the time step chosen for this splitting, as well as the curvature of 
the lines. 

We have demonstrated the flexibility of this method in four numerical examples. First 
we have shown that the time until a molecule binds to an operator site on DNA depends 
on the length between two road blocks on each side of the operator site. As the length 
increases, the time until a molecule binds decreases. The decrease levels off as the distance 
gets bigger than the average sliding distance of a molecule on the DNA, in agreement with 
Refs. I2S1 and ESI In the second example we consider two spirals embedded in a sphere. 
Molecules in space undergo a reversible reaction with the spirals and molecules on the 
spirals also react reversibly. In a third example we have considered a simple model of 
active transport. Molecules in a sphere react reversibly with embedded lines and move 
deterministically towards the center, while bound. The lines also move in space and we have 
shown how the radial distribution of the molecules is shifted due to the active transport. In 
a final example we have simulated lines that are growing and shrinking in space and shown 
how the dynamics of the lines correlates to the concentration of the molecules in the system. 
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